module bcs_mod
    !> 引入双精度常量定义
    use iso_fortran_env, only: real64
    use iso_fortran_env, only: int64
    implicit none
contains
    !> 使用划行划列法计算加入边界条件后的整体刚度矩阵 K' 和节点力向量 F'
    subroutine add_bcs_by_reorder(gk, bcs, forces, k_, f_)
        !> 参数说明:
        !>   gk       : 输入的整体刚度矩阵 (n×n)，奇异矩阵
        !>   bcs     : 边界条件数组 (m×2)，每行代表一个边界条件
        !>             第1列: 自由度索引 (1-based)，第2列: 边界条件值
        !>   forces  : 节点力向量 (p×2)，每行代表一个节点力
        !>             第1列: 自由度索引 (1-based)，第2列: 力值
        !>   k_      : 输出的刚度矩阵 (n×n)，非奇异
        !>   f_      : 输出的节点力向量 (n)
        !
        ! 注意:
        !   - 所有数组维度需在调用前正确分配
        !   - 使用单精度 (real) 计算
        real(real64), intent(in)  :: gk(:, :), bcs(:, :), forces(:, :)
        real(real64), intent(out) :: k_(:, :), f_(:)
        real(real64)    :: Uval
        integer :: i, j, k, ndof  !! 作为数组索引,不能够用int64

        ! 初始化
        ndof = size(gk, 1)  ! 矩阵维度
        k_ = gk             ! 复制刚度矩阵
        f_ = 0.0_real64            ! 初始化力向量
        print *, '  add_bcs_by_reorder->copy gk ok !'
        ! 处理节点力 (1-based 索引)
        do i = 1, size(forces, 1)
            j = int(forces(i, 1))  ! 直接使用 1-based 索引
            f_(j) = f_(j) + forces(i, 2)
        end do
        print *, '  add_bcs_by_reorder->add global forces to f_ array ok !'

        ! 处理边界条件 (1-based 索引)
        print *, '  ndof = ',ndof
        print *, '  k_ size=',shape(k_)
        print *, '  f_ size=',shape(f_)
        do i = 1, size(bcs, 1)
            j = int(bcs(i, 1))     ! 直接使用 1-based 索引
            Uval = bcs(i, 2)
            if (j>size(k_,1) .or. j>size(f_,1)) print *,'==> error index j over range:: add_bcs_by_reorder, j= ',j,' <=='
            
            do k = 1, ndof
                f_(k) = f_(k) - k_(k, j) * Uval
                if (k>size(k_,1)) then
                    print *,'==> error index k over range:: add_bcs_by_reorder, k= ',k,' <=='
                end if
            end do
            k_(:, j) = 0.0_real64  ! 整列清零 (Fortran 数组切片)
            k_(j, :) = 0.0_real64  ! 整行清零
            k_(j, j) = 1.0_real64  ! 对角线设为1
            f_(j) = Uval
        end do
        print *, '  add_bcs_by_reorder->add bcs to k_ and f_ array ok !'
    end subroutine add_bcs_by_reorder

    !> 使用划行划列法计算加入边界条件后的整体刚度矩阵 K' 和节点力向量 F'
    subroutine add_bcs_by_bignumber(gk, bcs, forces, k_, f_)
        !> 参数说明:
        !>   gk       : 输入的整体刚度矩阵 (n×n)，奇异矩阵
        !>   bcs     : 边界条件数组 (m×2)，每行代表一个边界条件
        !>             第1列: 自由度索引 (1-based)，第2列: 边界条件值
        !>   forces  : 节点力向量 (p×2)，每行代表一个节点力
        !>             第1列: 自由度索引 (1-based)，第2列: 力值
        !>   k_      : 输出的刚度矩阵 (n×n)，非奇异
        !>   f_      : 输出的节点力向量 (n)
        !
        ! 注意:
        !   - 所有数组维度需在调用前正确分配
        !   - 使用单精度 (real) 计算
        real(real64), intent(in)  :: gk(:, :), bcs(:, :), forces(:, :)
        real(real64), intent(out) :: k_(:, :), f_(:)
        integer :: i, j, ndof
        real(real64)    :: Uval,c

        ! 初始化
        c=10000000000.0_real64
        ndof = size(gk, 1)  ! 矩阵维度
        k_ = gk             ! 复制刚度矩阵
        f_ = 0.0_real64           ! 初始化力向量
        ! 处理节点力 (1-based 索引)
        do i = 1, size(forces, 1)
            j = int(forces(i, 1))  ! 直接使用 1-based 索引
            f_(j) = f_(j) + forces(i, 2)
        end do
        ! 处理边界条件 (1-based 索引)
        do i = 1, size(bcs, 1)
            j = int(bcs(i, 1))     ! 直接使用 1-based 索引
            Uval = bcs(i, 2)

            k_(j, j) = c+k_(j, j) ! 对角线元素乘大数
            f_(j) = f_(j) - c*Uval
        end do
    end subroutine add_bcs_by_bignumber
end module bcs_mod

